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ABSTRACT 


This  report  describes  in  detail  three  special  techniques  which 
can  be  specifically  applied  in  the  SPACETMCK  system  to  deter¬ 
mine  the  motion  of  an  artificial  satellite  under  the  influence 
of  perturbing  accelerations.  These  are  Variation  of  Parameters, 
Crowell's  Method,  and  Encke's  Method. 
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1.  INTRODUCTION 


The  purpose  of  this  report  is  to  describe  in  some  detail  three  special 
perturbation  techniques  vhich  find  specific  application  in  the  SPACETRACK 
system  for  determining  the  motion  of  an  artificial  satellite  under  the 
influence  of  perturbing  accelerations.  The  methods^  and  a  brief  comment  on 
each,  are: 

Variation  of  Parameters^  also  commonly  called  the  variation  of  elements , 
expresses  the  perturbations  in  the  orbit  as  a  function  of  the  elements, 
that  is  to  say,  the  difference  between  the  elements  of  the  orbit  at  epoch 
and  those  of  the  osculating  orbit  at  any  time  (t).  Variation  of  parameters 
is  used  in  SPWDC  and  SPIRDEC. 

Cowell ’ s  method  integrates  the  equations  of  motion  in  rectangular  coordinates 
directly,  giving  the  rectangular  coordinates  of  the  perturbed  body.  This 
method  finds  application  in  ESPOD  and  the  final  phase  of  SPIRDEC. 

Encke  * s  method  differs  from  Cowell's  in  that  it  is  the  difference  between 
where  the  body  actually  is  in  rectangular  coordinates  at  some  instant  and 
where  it  would  have  been  if  no  perturbative  accelerations  were  present 
(two  body  reference  orbit)  that  is  calculated.  This  method  is  used  in 
MUNENDC. 

The  remainder  of  this  report  is  devoted  to  a  development  of  these  three  methods. 
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2.  VARIATION  OF  PARAMETERS  METHOD 


a .  Description 

In  Figure  let  A  be  the  position  of  the  satellite  at  some  time  t^.  Suppose 

at  this  instant  all  the  perturbing  forces  acting  on  it  ceased  to  exist.  The 

satellite  would  then  continue  to  move 

in  an  elliptical  orbit  (AB)  with  the 

constant  elements  of  the 

classical  two  body  problem.  Further^ 

if  the  position  and  velocity  of  the 

satellite  is  known  at  t^^  we  have  the  requisite  number  of  equations  for 

determining  the  six  elements  of  the  ellipse.  The  ellipse  AB  is  called  the 

osculating  ellipse  at  time  t^  and  the  constants  C^j  *  * 

osculating  elements  at  time  t^.  The  actual  path  of  the  satellite  when  the 

effect  of  the  perturbations  is  taken  into  account  is  not  along  the  elliptical 

path  AB  but  rather  along  some  other  path  represented  by  AC  in  Figure  1.  The 

satellite  at  the  instant  t  has  the  same  coordinates  and^  by  definition, 

o 

the  same  velocity  components  in  the  unperturbed  as  in  the  perturbed  orbit. 
Stated  another  way,  the  satellite  has  the  position  and  is  moving  instant¬ 
aneously  as  it  would  in  purely  two-body  motion.  Obviously  one  could  com¬ 
pute  a  set  of  elements  to  define  an  osculating  orbit  at  any  point  of  the 

actual  orbit.  At  time  t.. ,  just  subsequent  to  t  ,  there  could  be  defined 

11  1  °  . 

a  new  set  of  osculating  elements  C^,  ’’’^6^  associated  with  the  cor¬ 

responding  position  and  velocity  af  the  satellite  at  point  C  in  its  actual 
path.  The  satellite’s  position  would  have  been  at  B  at  time  t^,  in  the  absence 
of  all  perturbations.  It  is  these  differences  between  the  elements, 

-  C^,  etc.  that  are  the  perturbations  of  the  elements  in  the 

interval  t^  -  t  . 

1  o 


*  It  should  be  noted  that  C^,  are  identical  in  the  Keplerian  case 

to  a,  e,  i,  uu.  Cl  and  T  or  some  combination  of  these  elements. 
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The  major  perturbation  encountered  in  dealing  with  the  motion  of  artificial 
satellites  is  the  acceleration  caused  by  the  oblateness  of  the  earth  which 
is  much  greater  than  the  perturbing  accelerations  produced  on  it  by  other 
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The  effect  is  that  the  elements  of  the 


bodies  in  the  solar  system. 


Keplerian  orbit  vary  as  a  function  of  time.  Here  the  line  of  nodes  and  the 
perigee  point  move  very  rapidly  under  the  noncentral  force  field.  When 
the  rates  of  change  of  the  elements  are  known,  the  future  orbital  characteristics 
of  the  satellite  can  be  predicted. 

b .  Analysis 

The  analysis  known  as  the  variation  of  parameters  begins  with  the  assumption 
that  the  cartesian  coordinates  defining  the  position  of  the  satellite  are 
known.  In  vector  notation  these  are 


(1) 


where  r=xi+yj+zk  and  j  and  k  are  unit  vectors  along  the 
Y  and  Z  axes  respectively.  The  equations  of  motion  of  a  satellite 
of  mass  M  under  the  central  attraction  of  the  earth  (mass  M^)  and  acted 
upon  by  a  disturbing  function  R  can  be  written  as 


(2) 


SR  i^3R  i9Rk  ,  /XM  ^  XM  \ 

where  V  R  =  ^  ^  and  u  =  k  (M  +  M  ). 

O  X  O  I  o  ^  o 


(3) 


From  equation  (l)  with  the  C  's  (k  =  1,  2,  "^...,6)  as  functions  of  time 

6  ^ 


(4) 


1.  Additional  perturbations  associated  with  earth  satellites  arise  from 
atmospheric  drag,  solar  radiation,  etc. 
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But  B  r  _  d  r  which  implies 
d  t  ■  d  t 


d  r  C 
S  C, 


k  =  0 


(5) 


k=l 


Invoking  (5)  and  differentiating  equation  (i<-)  with  respect  to  t,  we  obtain 

6 

(6) 


—  dr 


r  = 


4.2  dt  d  C, 

o  t  k=l  k 


Substituting  this  result  back  into  equation  (2)  yields 

^  ^k  =  V  R 

d  t^  r^  k=l  ^  *"k 


(7) 


For  the  osculating  orbit,  V  R  =  0  and  the  C  's  are  constants  so  that 


::>2  - 

+  ILJL  _  0 
^2  3  ■ 
St  r 

hence  from  equation  (t) 

A 


(8) 


2 


d^  r  C. 


k=l  ^^k 


k  =  V  R. 


(9) 


It  is  common  to  rewrite  equation  {9)j  making  use  of  the  fact  that 

d  I  d  r  ]  dr 


d^  r  _  _ 

dt  dC,  ~  3C,  \dt  /  Sc  ^ 
k  k  '  k 


as 


d  r  C 
S  C, 


k  =  V  r; 


(10) 


The  time  derivatives  of  the  orbital  elements  can  now  be  found  by  solving 
equations  (5)  and  (lO)  simultaneously  for  the  However,  the  solution 

of  these  equations  can  be  performed  more  readily  by  a  rearrangement  which 
introduced  new  functions  of  the  parameters  called  Lagrangian  brackets. 
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d  r 


dot  product  of  equation  (lO)  vith  equation  (5)  with 


If  we  take  the 
S  r 

and  subtract  the  tvo^  the  resulting  six  equations  may  be  written  as 

Cj 


k=l 


a  c.  'he, 

_  J  k 


h  r 

a  c. 


a  r 

a  c. 


^k  =  V  R  . 


a  r  (j=l,  2,.. .,6).  (ll) 


a  c. 

0 


The  quantity  in  brackets  in  equation  (ll)  is  Lagrange's  bracket  and  is 
commonly  denoted  by 


[fj,  ''3'  5  S 


where 


II 


5  X  5  X 

a  c.  a  c, 


jj;  ^  (z;  il  ^ 


(la) 


a  X 


a  X 


a  c.  a  c, 

J  k 


with  similar  expressions  for  a  a~'[c^  ^  )  ‘ 

^  j,  k-*  ^  y 

The  right  hand  side  of  equation  (ll)  is  the  partial  derivative  of  R  with 

respect  to  C  .  which  is 
J 


a  R  a  X  .  a  R  ^  ^  ^  ^  =  a  r  (13) 

ax  a  c.  §  a  c .  a  z  a  c.  a  c . 

j  y  0  j  j 

Using  equations  (12)  and  (l3)^  equation  (ll)  may  be  written  very  simply  as 


a  R  (j=l,  2,. ..,6). 
a  c. 

J 


(1^+) 


It  is  these  six  equations  that  are  to  be  solved  for  C^.  Equation  (l^) 
contains  thirty-six  Lagrangian  brackets  but  from  the  definition  of  these 
brackets  in  equation  (l2)  it  is  noted  that 


5,  ^2  -  °  Ej.  ■  '&■  "2  ■ 

Therefore,  the  number  of  distinct  Lagrangian  brackets  to  be  evaluated  is 
only  fifteen  instead  of  the  original  thirty-six. 
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An  example  of  the  differential  equations  representing  the  variation  of 

orbital  elements  which  results  from  evaluating  the  Lagrangian  brackets 

.  (2) 

are: 


da 

-  _2_  SR 

dt 

na  5M' 

de 

1  -  SR 

/l  -  e‘ 

■  SR 

dt 

2  dM 

2 

Siu 

na  e 

na  e 

dou 

cos  i 

dt 

2  T  ’ 

na  y/1  - 

e  sin  1 

Si 

- 

cos  i 

SR 

dt 

2  f.  5 

na  yi  -  e 

sin  i 

Sou 

^  - 

1 

dt 

2  y, 

na  y  1  -  e 

sin  'i 

Si 

m  _ 

1  - 

dR _ 2 

SR 

dt 

n-  2 

na  e 

he  na 

3a 

2 

na  e 


he 


The  equations  above  differ  from  the  differential  equations  solved  for  in 
SPWDC  and  SPIRDEC  vhich  are  in  terms  of  an  N-M  element  set.  In  terms  of  the 
N-M  element  set  and  considering  perturbations  due  to  the  earth’s  bulge, 

(3) 

radiation-pressure  and  drag,  the  variation  of  parameter  equations  become 


(2)  Kozai,  Y.,  ”The  Motion  of  a  Close  Earth  Satellite",  The  Astronomical 
Journal,  64,  No.  9^  Page  3^9;  November  1959* 

(3)  Aeronutronic,  Special  Perturbations  Weighted  Differential  Correction 
Program  Document,  Technical  Documentary  Report  No.  ESD-TDR- 63-645^  11  Dec.  63* 
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In  SPWDC, 
technique 


dt 


k  L  n 
e  + 


) 


e  — 


) 


—  k  h 
dt  “  ® 


N 


these  equations  are  solved  using  a  Runge-Kutta  numeric  integration 
to  determine  the  satellite's  position  and  velocity. 
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3.  COWELL  *S  METHOD 


a.  Description 

Covell’s  method  of  numerical  integration  makes  no  explicit  use  of  a  conic 
section  as  the  first  approximation  to  the  orbit^  but  rather^  the  equations 
of  motion  in  rectangular  coordinates  are  integrated  directly^  giving  the 
rectangular  coordinates  of  the  disturbed  body.  The  origin  is  usually  taken 
at  the  primary  body^  but  this  restriction  is  not  necessary  since  the  center 
of  mass  of  the  system  or  of  any  of  the  disturbing  bodies  may  be  used.  The 
only  restriction  is  that  the  motion  of  all  bodies  exerting  appreciable 
effects  are  knovn  relative  to  the  chosen  origin  at  some  time.  The  only 
practical  disadvantage  of  the  method  is  that  the  integrals  contain  many 
significant  figures  and  change  rapidly  vith  time.  In  consequence,  the 
integration  tables  are  slovly  convergent  vhich  compels  the  use  of  a  small 
tabular  interval. 

b .  Analysis 

Consider  tvo  point  masses,  m  and  m^, 
with  coordinates  and 

relative  to  a  Cartesian  coordinate 
system  X,  Y  and  Z  as  illustrated  in 
Figure  2 . 


Let  r  be  a  vector  from  the  origin  of 
a 

the  coordinate  system  to  m  and  r^  be 

a  b 


a  vector  to  m^. 

In  addition,  let  s  be  a  unit  vector  in  the  direction  and  i,  j,  and  k 

unit  vectors  from  the  origin  along  X,  Y  and  Z. 


From  Nevton*s  lav,  the  gravitational  attraction  betveen  m^  and  m^  can  be 
expressed  by: 


F  = 


,2 

k  m^  m^ 


(l6) 
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where  r  is  the  distance  between  and  and  k  is  a  constant  of  proportionality 
depending  on  the  units  of  mass,  time  and  distance  chosen.  The  force  on  m 


due  to 


IS 


k  m 


F  =  m  r  = 
a  a  a 


a  \  " 

2 -  ^ 


(it) 


and  similarly,  that  on  m^  due  to  m^  is 


"'b  ”  "b  ''b 


k  m^  m^  s 


(18) 


To  determine  the  components  of  the  force  acting  on  m  in  the  X,  Y,  Z 

a 

directions,  i.e.,  the  r|-,  and  Q-  components,  it  is  necessary  to  obtain 
the  dot  product  of  F  with  the  unit  vectors  i,  j  and  k.  Thus,  the  com¬ 

ponent  of  the  force  on  m  is 

a 


-  -  -'2  ^a^ 

F  ^  F  .i  =  m  r  .i  =  m  E!  =  k  — r —  cos  (F  ,  i),  or, 
a?  a  a  a  a  ^a  2  ^  a^  ^  ^ 


m 


2 

?  =  k  m  m, - r - 

^a  a  b  3 


where  r  = 


(?a  -  5b "  (^a  -  *  <5a  ’  Cb^  ' 

Similarly,  the  5“  component  of  the  force  on  ra^  due  to 


m  IS 


"b  5b  -  "  “b 


(5a  -  Eb) 


m 


(19) 


(20) 


Similar  expressions  can  be  written  for  the  r|-  and  components. 


If  additional  point  masses  m^,  m^,  m^. 


are  introduced  into  this  system 


and  denoting  any  one  of  these  point  masses  by  m.,  equations  similar  to  (l9) 

J 

and  (20)  expressing  the  total  accelerations  of  m^  and  m^  may  be  obtained  by 
summing  all  these  attractions. 
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In  figure  2^  represents  one  such  mass  vhose  distance  from  m^  and  m^^  is 

l,b 


and  Pt  respectively 


The  expression  then  for  the  §-  components  would  be 

2  ^^b  ■  ^ 

m  ?  =  k  m  m,  - = — 

a  ^a  an  3 


(21) 


-  2  ■  ^b'  ^ 

”b  ^b  =  ^  “b”a  - 3 -  ^  4 


,2  ■  ^b^ 

^"b”j  '73 - 


J,  b 


vhere 


J,a 


J,h 


(^s  -  ^  \  -  V'  ^  -  Cj 


g 


<h  -  *  K  -  ■ij’'  *  <‘b  -  Cj)^  ^ 


(22) 

(23) 

(2k) 


Again,  similar  expressions  can  be  written  for  the  r)-  and  Q-  components.  Let 
the  origin  of  coordinates  be  taken  at  m  which  is  equivalent  to  the  linear 

Q 

transformation 

=  x;  5  .  -t  =  X .  • 

It  then  follovs  directly  from  (25)  that 


(25) 


-5b  =  ■ 


(26) 


Finally,  put 


2  2  ^  2  ^  2 

r  =  X  +  V  +  z 
j  j  j 


;Pj=  |U 


J  *  <"j 


(2T) 
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Divide  equation  (2l)  by  and  equation  (22)  by  and  then  subtract  (2l)  from 
(22).  The  result  is  the  equation  of  motion  of  m^^  relative  to  m^  and  can  be 
expressed  as 


X  =  -  k  (m  + 
a 


^-^,2  “j  ^ "  Z- 

J 


2 

J 


(28) 


A  more  familiar  expression  for  equation  (28)  is 


X  =  -  k 


X  Y”  p  ^ 

(“a  "  "  Z_  "  "j  4"  ■; 

J  \ 


(29) 


This  equation  and  similar  equations  in  y  and  z  are  the  fundamental  equations 
in  Cowell's  method.  Other  accelerations  can  be  combined  in  the  right  hand 
side  of  equation  (29)  such  as  drag^  zonal  harmonics  of  the  earthy  radiation 
pressure^  etc. 
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4.  ENCKE'S  METHOD 


a.  Description 

Encke^s  method  differs  from  Coveil's  in  that  the  coordinates  of  the  disturbed 
body  are  not  obtained  directly  but  rather  from  the  difference  betveen  the 
position  the  body  vould  have  in  an  osculating  orbit  referenced  to  some  time  t, 
called  the  epoch  of  osculation^  and  its  true  position  that  is  calculated. 

The  departures  from  the  osculating  orbit  are  the  perturbations.  The  advantage 
of  this  method  is  that  for  times  near  the  epoch  of  osculation^  the  perturbations 
are  small  and  can  be  expressed  by  a  fev  significant  figures  permitting  a 
larger  tabular  interval  than  vith  Covell's  method.  Against  this  is  the 
disadvantage  that  each  step  of  Encke's  method  takes  longer  and  the  perturbations 
increase  vith  time  requiring  an  occasional  redetermination  of  the  osculating 
orbit. 

b .  Analysis 

Let  m^  and  nij^  be  tvo  point  masses  vith  y^  z  the  coordinates  of  m^  relative 
to  m^  vhere  is  moving  under  the  attraction  of  m^  alone.  The  equations  of 
motion  are  knovn  to  be 


(30) 


(31) 


(32) 


vhere 


r 


o 
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CV;  P;  Y  represent  the  perturbations  produced  by  the  presence  of  additional 
masses  such  that  the  true  position  (x^  and  z)  of  at  any  time  can  be 
represented  by 


X  =  +  a,  y  =  y^  +  p,  z  =  +  Y, 


and  the  actual  equations  of  motion  are 
2 


X  = 


m 


J 


J 

..  1. 


" 


X  . 


(33) 


(34) 


with  similar  expressions  for  y  and  *z\  Note  that 

/2  2  ^  \  ^ 
r  =  (x  +  y  +  z;*^. 


Subtracting  (30)  from  (3^)  yields 


2 

x-x^  =a’  =  k  (nij^  +  m^) 


\r, 


m 


J 


J 


fx.  -  X  X, 

_J _  - 

3  3 

\Pj  "J  / 


(35) 


again  with  similar  expressions  for  p  and  Y- 

Equation  (35)  could  be  solved  for  a,  P,  y  direct  integration  by  calculating 

X  X 

o  for  each  step  by  the  laws  of  elliptic  motion  and  the  term  —  at  each  step 

r  3 
o 

by  extrapolating  a  and  adding  it  to  x  to  give  x,  etc.,  but  this  approach  would 
not  be  convenient  in  practice  for,  since  a  is  a  small  quantity,  o  is  nearly 

r  3 

X  o 

equal  — ,  and  these  two  terms  would  have  to  be  calculated  to  many  more 

significant  figures  than  are  needed  in  their  difference.  Encke  developed 
the  following  transformation  to  overcome  this  difficulty.  Treating  only 
the  equation  for  5,  since  those  for  ^  and  y  exactly  similar. 


r3  r 
o 


X  1 

“3 "  7 


O  - 


xr  I  1 
o  I  = 

1-  - 

/  r- 


X  -  a 


(36) 


1.  Refer  to  equation  (29) 
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But 


=  (x^  +  a)^  +  (y^  +  +  (z^  +  y)^ 

=  r  ^  +  2x  a  +  2y  P  +  2z  Y  ■'■  3^  ■*■ 

o  o  o  o  ^  ^ 


Dividing  ( 37 )  by  r  yields 


(37) 


1  +  2 


(Xq  ■*■  i  Q')a+  (y^  +  ^)P  +  (z^  +  |y)y 


and  putting 


q  = 


(Xq  |a)Qf  +  (Yq  +  ?3)3  +  (Zq  +^y)y  ^ 


(38) 


(39 


equation  (38)  may  be  rewritten  as 


=  1  +  2  q. 


('^O) 


From  this^ 


V  =  (1  - 


=  1  -  (1  +  2  q) 


r3A 


(i^l) 

(ii2) 


If  we  now  define  a  function  f  by 

^  1  -  (1  +2 
q 

equation  (35)  can  be  rewritten  in  the  form 


f  = 


a  =  ("S  +  "^a)  (f  q  X  -  a)  +  ^ 


1 2  /X. 

km./  ,1  -  X 

p  3 


r3 

J 


Similar  equations  for  ^  and  V  can  be  expressed;  it  is  these  equations  that 
are  solved  in  Encke*s  method.  Again,  as  in  the  discussion  of  Cowell’s  method, 
additional  perturbing  accelerations  can  be  added  to  the  right  half  of  equation 


lii 
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